Import libraries¶

Install them with pip install -r requirements.txt

In [ ]:
import urbanpy as up
import contextily as cx
import matplotlib.pyplot as plt
import plotly
import plotly.express as px
from mpl_toolkits.axes_grid1 import make_axes_locatable
from tqdm.auto import tqdm

/Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm
In [ ]:
# Activate tqdm progress bar for pandas.apply
tqdm.pandas()

Get city administrative limits¶

In [ ]:
pos_id = 0  # Result position
riyad = up.download.nominatim_osm("Riyadh Governorate, Saudi Arabia", pos_id)  # Nominatim query
In [ ]:
# Plot administrative limits
ax = riyad.plot(facecolor="none", edgecolor="r")
# Add a basemap
cx.add_basemap(ax, source=cx.providers.CartoDB.Positron, crs="EPSG:4326")
ax.set_axis_off()
No description has been provided for this image

Get city population density¶

In [ ]:
# Query the population data in Saudi Arabia from the repository of Meta Data For Good Population Maps in the Humanitarian Data Exchange platform
population_search_results = up.download.search_hdx_dataset("Saudi Arabia")
In [ ]:
population_search_results
Out[ ]:
created name population size_mb url
id
1 2020-03-05 sau_general_2020_csv.zip Overall population density 45.71 https://data.humdata.org/dataset/14b288ca-1855...
8 2019-09-23 sau_children_under_five_2020_csv.zip Children (ages 0-5) 45.55 https://data.humdata.org/dataset/14b288ca-1855...
9 2019-09-23 sau_elderly_60_plus_2020_csv.zip Elderly (ages 60+) 45.48 https://data.humdata.org/dataset/14b288ca-1855...
10 2019-09-23 sau_men_2020_csv.zip Men 45.66 https://data.humdata.org/dataset/14b288ca-1855...
11 2019-09-23 sau_women_2020_csv.zip Women 45.63 https://data.humdata.org/dataset/14b288ca-1855...
12 2019-09-23 sau_women_of_reproductive_age_15_49_2020_csv.zip Women of reproductive age (ages 15-49) 45.62 https://data.humdata.org/dataset/14b288ca-1855...
13 2019-09-23 sau_youth_15_24_2020_csv.zip Youth (ages 15-24) 45.52 https://data.humdata.org/dataset/14b288ca-1855...
In [ ]:
# We selected the index 9 that refers to the Elderly (ages 60+)
saudi_arabia_pop = up.download.get_hdx_dataset(population_search_results, 9)
saudi_arabia_pop.head()
Out[ ]:
longitude latitude sau_elderly_60_plus_2020
0 39.847083 33.000139 0.141545
1 40.316806 33.000139 0.141545
2 44.159583 33.000139 0.367002
3 44.166250 33.000139 0.178008
4 44.197083 33.000139 0.178008
In [ ]:
# Number of population points for the country
saudi_arabia_pop.shape[0]
Out[ ]:
9893484
In [ ]:
# Use the city adm limits to filter the country population data
filtered_pop = up.geom.filter_population(saudi_arabia_pop, riyad)
filtered_pop.head()
Out[ ]:
longitude latitude sau_elderly_60_plus_2020 geometry
4699315 46.770417 25.267083 0.035855 POINT (46.77042 25.26708)
4699316 46.805694 25.267083 0.035855 POINT (46.80569 25.26708)
4699317 46.805972 25.267083 0.035855 POINT (46.80597 25.26708)
4699318 46.806528 25.267083 0.035855 POINT (46.80653 25.26708)
4699319 46.806806 25.267083 0.035855 POINT (46.80681 25.26708)
In [ ]:
# Number of population points for the city
filtered_pop.shape[0]
Out[ ]:
529876
In [ ]:
# Plot population points
ax = filtered_pop.plot("sau_elderly_60_plus_2020", markersize=0.01, legend=True)
# Plot administrative limit
riyad.plot(facecolor="none", edgecolor="r", ax=ax)
# Add a basemap
cx.add_basemap(ax, source=cx.providers.CartoDB.Positron, crs="EPSG:4326")
ax.set_axis_off()

/Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/geopandas/plotting.py:732: FutureWarning: is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
  if pd.api.types.is_categorical_dtype(values.dtype):
No description has been provided for this image

To improve the interpretability of the population plot we will group them in uniform spatial units known as H3 Hexagons.¶

In [ ]:
# Generate hexagons of resolution 7 (~5.1612km2)
riyad_hexs = up.geom.gen_hexagons(resolution=7, city=riyad)

See the table of resolution and sizes here

In [ ]:
# Plot the city H3 hexagons
ax = riyad_hexs.plot(facecolor="none", edgecolor="r")
# Add a basemap
cx.add_basemap(ax, source=cx.providers.CartoDB.Positron, crs="EPSG:4326")
ax.set_axis_off()
No description has been provided for this image
In [ ]:
# Sum the point's population within each hexagons
merges_hexs = up.geom.merge_shape_hex(
    riyad_hexs, filtered_pop, {"sau_elderly_60_plus_2020": "sum"}
)  # You can use other aggregation methods like max, min, count, and mean
In [ ]:
# Plot hexagons colored by population density
ax = merges_hexs.plot("sau_elderly_60_plus_2020", legend=True, missing_kwds={"color": "grey"})
# Add a basemap
cx.add_basemap(ax, source=cx.providers.CartoDB.Positron, crs="EPSG:4326")
ax.set_axis_off()

/Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/geopandas/plotting.py:732: FutureWarning: is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
  if pd.api.types.is_categorical_dtype(values.dtype):
No description has been provided for this image
In [ ]:
# Get amenities (Point of Interest) from Overpass
gdf_pois, _ = up.download.overpass(type_of_data="node", query={"amenity": ["clinic", "hospital"]}, mask=riyad)
gdf_pois.head()
Out[ ]:
type id lat lon tags geometry poi_type
38 node 4378219310 24.400960 46.835173 {'amenity': 'hospital', 'name': 'مستوصف الحاير... POINT (46.83517 24.40096) hospital
7 node 1354660020 24.583313 46.864015 {'amenity': 'hospital'} POINT (46.86402 24.58331) hospital
73 node 6953249575 24.678829 46.776056 {'addr:city': 'الرياض', 'addr:housenumber': 'ح... POINT (46.77606 24.67883) hospital
3 node 746341536 24.687040 46.797248 {'amenity': 'hospital'} POINT (46.79725 24.68704) hospital
74 node 7011337045 24.694442 46.790222 {'addr:street': 'الزبير بن العوام', 'amenity':... POINT (46.79022 24.69444) hospital

See all types of data you can query from Overpass

In [ ]:
# Plot points of interests
ax = gdf_pois.plot("poi_type", legend=True, figsize=(10, 10))
# Add a basemap
cx.add_basemap(ax, source=cx.providers.CartoDB.Positron, crs="EPSG:4326")
ax.set_axis_off()

/Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/geopandas/plotting.py:732: FutureWarning: is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
  if pd.api.types.is_categorical_dtype(values.dtype):
No description has been provided for this image
In [ ]:
# Start routing server (needs docker)
up.routing.start_osrm_server("gcc-states", "asia", "car")  # foot,car,bicycle

Starting server ...
osrm_routing_server_asia_gcc-states_car
Server was started succesfully
In [ ]:
# Calculate travel times (duration in minutes and distance in km) from hexagons to points of interest
merges_hexs_tt = up.accessibility.travel_times(merges_hexs, gdf_pois)

/Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/urbanpy/accessibility/accessibility.py:260: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  gdf["lon"] = gdf.geometry.centroid.x
/Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/urbanpy/accessibility/accessibility.py:261: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  gdf["lat"] = gdf.geometry.centroid.y
  0%|          | 0/1627 [00:00<?, ?it/s]
Waiting for server to be ready ...

100%|██████████| 1627/1627 [00:15<00:00, 106.13it/s]
In [ ]:
up.routing.stop_osrm_server("gcc-states", "asia", "car")

Server was stoped succesfully
In [ ]:
merges_hexs_tt.head()
Out[ ]:
hex geometry sau_elderly_60_plus_2020 lon lat nearest_poi_ix distance_to_nearest_poi duration_to_nearest_poi duration_to_nearest_poi_label
0 875355da4ffffff POLYGON ((46.60462 24.61013, 46.60145 24.59749... 739.076001 46.614256 24.601532 51 8.2799 10.111667 0-15
1 8753734caffffff POLYGON ((46.50737 24.92900, 46.50419 24.91635... 51.435090 46.517054 24.920377 96 25.7935 22.543333 15-30
2 875355891ffffff POLYGON ((47.00107 24.54897, 46.99787 24.53637... 1.507293 47.010623 24.540370 1 22.5693 46.918333 45-60
3 875355809ffffff POLYGON ((47.21802 24.66309, 47.21479 24.65051... 7.536465 47.227535 24.654486 49 65.6733 111.828333 90-120
4 87535582effffff POLYGON ((47.17685 24.73132, 47.17362 24.71873... NaN 47.186380 24.722705 49 61.4836 101.770000 90-120
In [ ]:
# Plot distance and duration maps
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 10))

# Plot distance
divider = make_axes_locatable(ax1)
cax = divider.append_axes("right", size="5%", pad=-0.1)
merges_hexs_tt.plot("distance_to_nearest_poi", legend=True, ax=ax1, cax=cax)
cx.add_basemap(ax1, source=cx.providers.CartoDB.Positron, crs="EPSG:4326")
ax1.set_axis_off()
ax1.set_title("Distance to Nearest PoI")
# Plot duration
merges_hexs_tt.plot("duration_to_nearest_poi_label", cmap="magma_r", legend=True, ax=ax2)
# Add a basemap
cx.add_basemap(ax2, source=cx.providers.CartoDB.Positron, crs="EPSG:4326")
ax2.set_axis_off()
ax2.set_title("Duration to Nearest PoI")

plt.show()

/Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/geopandas/plotting.py:732: FutureWarning: is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
  if pd.api.types.is_categorical_dtype(values.dtype):
/Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/geopandas/plotting.py:732: FutureWarning: is_categorical_dtype is deprecated and will be removed in a future version. Use isinstance(dtype, CategoricalDtype) instead
  if pd.api.types.is_categorical_dtype(values.dtype):
No description has been provided for this image

Generate interactive maps¶

In [ ]:
plotly.offline.init_notebook_mode()
In [ ]:
fig = up.plotting.choropleth_map(
    merges_hexs_tt, "sau_elderly_60_plus_2020", title="Estimated Elderly (60+) Population - 2020", opacity=0.5
)

# Remove the hexagon outlines to make the map clearer
fig.update_traces(marker_line_width=0)

# Make space for the title
fig.update_layout(margin=dict(l=0, r=0, b=0))
In [ ]:
# Filter out the hexagons without population
merges_hexs_tt_filtered_pop = merges_hexs_tt.query("sau_elderly_60_plus_2020 > 0").reset_index(drop=True)
In [ ]:
# Get ordered category labels
category_orders = merges_hexs_tt_filtered_pop["duration_to_nearest_poi_label"].unique().sort_values()
In [ ]:
fig = up.plotting.choropleth_map(
    merges_hexs_tt_filtered_pop,
    color_column="duration_to_nearest_poi_label",
    color_discrete_sequence=px.colors.sequential.Plasma_r,
    category_orders={"duration_to_nearest_poi_label": category_orders},
    labels={"duration_to_nearest_poi_label": "Minutes"},
    title="Travel Time to Nearest PoI",
    opacity=0.5,
)

# Make space for the title
fig.update_layout(margin=dict(l=0, r=0, b=0))

# Remove the hexagon outlines to make the map clearer
fig.update_traces(marker_line_width=0)

/Users/claudio/Documents/urbanpy-riyad/.env/lib/python3.11/site-packages/plotly/express/_core.py:2044: FutureWarning:

The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning.
In [ ]: